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Abstract.  In  the  field  of  nondestructive  evaluation,  new  and  improved  techniques  are  constantly 
being  sought  to  facilitate  the  detection  of  hidden  corrosion  and  flaws  in  structures  such  as  airplanes 
and  pipelines.  In  this  paper,  we  explore  the  feasibility  of  detecting  such  damages  by  application  of 
an  eddy  current  based  technique  coupled  with  reduced  order  modeling. 

We  begin  by  developing  a  model  for  a  specific  eddy  current  method  in  which  we  make  some 
simplifying  assumptions  reducing  the  three-dimensional  problem  to  a  two-dimensional  problem.  (We 
do  this  for  proof-of-concept.)  Theoretical  results  are  then  presented  which  establish  the  existence  and 
uniqueness  of  solutions  as  well  as  continuous  dependence  of  the  solutions  on  the  parameters  which 
represent  the  damage.  We  further  discuss  theoretical  issues  concerning  the  least  squares  parameter 
estimation  problem  used  in  identifying  the  geometry  of  the  damage. 

To  solve  the  identification  problem,  an  optimization  algorithm  is  employed  which  requires  solving 
the  forward  problem  numerous  times.  To  implement  these  methods  in  a  practical  setting,  the  forward 
algorithm  must  be  solved  with  extremely  fast  and  accurate  solution  methods.  In  constructing 
these  computational  methods,  we  employ  reduced  order  Proper  Orthogonal  Decomposition  (POD) 
techniques.  This  approach  permits  one  to  create  a  set  of  basis  elements  spanning  a  data  set  consisting 
of  either  numerical  simulations  or  experimental  data.  We  discuss  two  different  algorithms  for  forming 
the  POD  approximations,  a  POD/Galerkin  technique  and  a  POD/Interpolation  technique. 

Finally,  results  of  the  inverse  problem  associated  with  damage  detection  are  given  using 
both  simulated  data  with  relative  noise  added  as  well  as  experimental  data  obtained  using  a 
giant  magnetoresistive  (GMR)  sensor.  The  experimental  results  are  based  on  successfully  using 
experimental  data  to  form  the  POD  basis  elements  (instead  of  numerical  simulations),  thus 
illustrating  the  effectiveness  of  this  method  on  a  wide  range  of  applications.  In  both  instances 
the  methods  are  found  to  be  efficient  and  robust.  Moreover,  the  methods  were  fast;  our  findings 
demonstrate  a  significant  reduction  in  computational  time. 


1.  Introduction 

Nondestructive  evaluation  (NDE)  is  the  process  of  examining  a  material  or  article  without 
impairing  its  future  usefulness.  NDE  is  sometimes  referred  to  as  nondestructive  testing  (NDT) 
or  nondestructive  inspection  (NDI),  although  there  may  be  subtle  differences  in  their  definitions 
depending  on  the  author.  For  the  purposes  of  this  paper,  however,  we  will  use  the  terminology 
interchangeably. 

The  process  of  examining  a  material  using  nondestructive  evaluation  techniques  is  not  new  but 
is  becoming  increasingly  important  as  technology  continually  advances.  According  to  the  American 
Society  of  Nondestructive  Testing,  the  term  NDT  includes  many  methods  that  can;  (i)  detect 
internal  or  external  imperfections,  (ii)  determine  structure,  composition,  or  material  properties,  or 
(iii)  measure  geometric  characteristics.  Some  typical  structures  or  products  inspected  through  the 
use  of  NDE  technology  are  components  of  airplanes,  motor  vehicles,  pipelines,  bridges,  trains,  and 
power  stations. 

Nondestructive  evaluation  techniques  can  be  broken  down  into  seven  main  categories:  (i) 
visual  inspection,  (ii)  liquid  penetration  inspection,  (iii)  radiography  and  radiation  testing,  (iv) 
electromagnetic  testing,  (v)  acoustic  emission  monitoring,  (vi)  magnetic  methods,  and  (vii) 
ultrasonic  testing.  Each  of  these  categories  is  broad,  containing  within  them  several  specific  testing 
techniques.  The  choice  of  an  appropriate  NDE  technique  depends  on  the  specific  application.  We 
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will  not  go  into  detail  on  the  benefits  and  limitations  of  each  method,  but  will  instead  refer  the 
reader  to  [8,  9,  18,  20,  38,  40]  for  a  thorough  discussion.  In  this  paper,  we  will  focus  only  on  a 
particular  electromagnetic  testing  method,  called  the  eddy  current  method. 

Eddy  currents  are  currents  found  in  any  conducting  material  which  is  subjected  to  a  time- 
varying  magnetic  field.  They  are  useful  for  NDE  purposes,  because  if  a  flaw  is  present  within  a 
conducting  material,  the  flow  of  the  eddy  currents  will  be  disrupted  in  some  manner.  From  this 
disruption,  and  its  effect  on  the  magnetic  flux  density,  one  can  often  discern  information  about  the 
damage  or  defect  within  the  material.  A  thorough  analysis  of  eddy  currents  and  their  behavior  can 
be  found  in  [44], 

Since  eddy  currents  are  only  found  in  a  conducting  material,  the  use  of  eddy  current  methods  for 
interrogation  is  limited;  however,  these  methods  have  been  proven  extremely  useful  when  examining 
pipeline  structures  and  aging  aircrafts.  There  are  many  devices  and  eddy  current  techniques  in  use 
today  including  the  self-nulling  eddy  current  probe  [49]  along  with  conformal  mapping  techniques 
[50],  the  magneto-optic/eddy  current  imager  [16,  45]  in  conjunction  with  eddy  current  imaging 
[17,  19],  the  SQUID  (Superconducting  Quantum  Interference  Device)  through  the  use  of  either 
injected  current  methods  or  induced  eddy  current  methods  [11,  13,  22,  37,  39,  48]  and  the  GMR 
(Giant  Magnetoresistive)  sensor  based  on  the  self- nulling  probe  design  [51,  52],  Each  of  these 
instruments  have  unique  features,  making  some  instruments  easier  to  use  than  others  or  more 
practical  depending  on  the  circumstances.  Some  of  the  instruments  provide  images  of  the  damage, 
while  others  provide  quantitative  data.  In  this  paper,  we  investigate  computational  methods  for  the 
estimation  of  a  damage  or  flaw  within  a  material  using  quantitative  data  taken  by  an  appropriate 
instrument.  In  some  cases,  instruments  providing  images  of  the  damage  may  give  us  valuable  a 
priori  information  about  the  defect,  allowing  us  to  obtain  a  more  accurate  estimate  of  the  damage. 

Damage  detection  is  quite  naturally  formulated  in  the  context  of  inverse  problems;  in  the  case 
of  electromagnetic  probes  these  involve  Maxwell’s  equations  at  some  level.  Consequently  one  can 
anticipate  computational  algorithms  that  are  time  and  computer  memory  intensive.  Motivated 
by  goals  of  on-line,  real-time  algorithms  to  be  used  in  portable  testing  devices,  one  is  led  to 
reduced  order  computational  ideas.  In  an  earlier  paper  [4]  (see  also  [3]),  we  suggested  techniques 
based  on  Proper  Orthogonal  Decomposition  (also  called  Principal  Component  Analysis)  methods 
and  gave  some  initial  computational  results  based  on  numerical  simulations  for  one-dimensional 
damages  which  suggested  significant  potential  for  such  an  approach.  In  this  paper  we  continue  our 
investigations,  providing  theoretical  foundations  for  the  proposed  approach  along  with  a  summary 
of  extensive  tests  using  numerical  simulations  for  one-dimensional  and  two-dimensional  damages. 
Moreover,  we  provide  a  summary  of  our  efforts  and  results  in  using  the  reduced  order  algorithms 
with  data  from  subsurface  damage  experiments  designed  and  carried  out  specifically  to  test  the 
efficacy  of  our  proposed  approach. 

We  first  analyze  a  specific  implementation  of  the  eddy  current  method  and  develop  a  model 
describing  the  behavior  of  the  magnetic  flux  density  as  a  relationship  to  the  total  current.  As  a 
consequence,  we  will  be  able  to  obtain  information  about  the  damage  from  the  relationship  between 
the  disruption  in  the  eddy  current  and  the  resulting  magnetic  flux  density. 
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As  discussed  in  the  previous  section,  eddy  currents  are  currents  found  in  any  conducting  material 
subjected  to  some  time-varying  magnetic  field  and  can  be  induced  in  a  variety  of  ways.  In  this 
paper,  we  limit  our  discussion  to  only  one  implementation:  we  examine  the  process  of  inducing 
eddy  currents  within  a  sample  by  placing  a  thin  conducting  sheet  carrying  a  uniform  current  above 
the  sample.  The  current  within  the  sheet  produces  a  magnetic  field  perpendicular  to  it  that  in 
turn  produces  eddy  currents  within  the  sample.  The  presence  of  a  flaw  within  the  sample  causes  a 
disruption  in  the  flow  of  the  eddy  currents  and  this  disruption  manifests  itself  in  the  magnetic  flux 
density  which  can  be  measured  by  a  device  placed  above  the  conducting  sheet.  A  schematic  of  the 
inspection  process  is  shown  in  Figure  1. 

By  assuming  uniformity  in  the  direction  of  the  current  flow  in  the  conducting  sheet  (labeled  the 
negative  z  direction,  denoting  the  coordinate  for  the  width  of  the  sample  in  Figure  2),  we  are  able  to 
reduce  the  three-dimensional  setup  described  above  to  a  two-dimensional  problem  in  the  xy  plane, 
where  x  denotes  the  coordinate  of  the  length  of  the  sample  and  y  denotes  the  coordinate  of  the 
thickness  of  the  sample.  We  make  this  simplifying  assumption  for  proof-of-concept,  to  illustrate  the 
feasibility  of  reconstructing  the  geometry  of  a  damage.  In  addition,  in  order  to  disregard  boundary 
effects  along  the  edges  of  the  sample,  we  assume  the  sample  and  conducting  sheet  are  infinitely  long, 
i.e.,  they  are  infinite  in  the  x  direction.  If  the  conducting  sheet  and  sample  are  not  of  infinite  extent, 
we  have  to  take  into  account  the  discontinuities  in  the  current  flow  at  the  boundaries.  Furthermore, 
the  damage  (which  we  shall  refer  to  as  a  “crack”)  is  assumed  to  be  rectangular  in  shape  and  centered 
along  the  length  of  the  sample  (along  the  x  direction)  at  x  =  0. 


Id 


Pick-up  coil  of  sensor 


2 


X1' 


Conducting  Sheet 


Sam  pie 


Material  Flaw 


Figure  1.  3-D  Schematic  of  Eddy  Current  Inspection  Process 

For  computational  purposes,  we  examine  a  finite  “window”  of  the  overall  problem,  which  is 
called  the  computational  domain  fh  In  choosing  the  boundaries  of  the  computational  domain  along 
the  length  of  the  sample  (x  boundaries),  we  recall  that  the  sample  and  conducting  sheet  are  assumed 
to  be  of  infinite  extent.  Therefore,  we  can  arbitrarily  choose  these  boundaries  by  assigning  evenly 
symmetric  boundary  conditions  to  account  for  the  infinite  extent  of  the  materials.  However,  since 
the  damage  is  centered  along  the  length  of  the  sample,  we  need  to  only  consider  half  of  the  sample 
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for  computational  purposes.  The  other  half  will  behave  in  the  same  manner.  Therefore,  we  choose 
the  left  boundary  to  be  located  at  the  center  of  the  crack  in  the  x  direction,  labeled  x  =  0,  with  the 
crack  symmetric  through  the  yz  plane  at  x  =  0,  and  the  right  boundary  is  chosen  at  x  =  50mm.  The 
y  boundaries,  or  top  and  bottom  boundaries,  are  at  y  =  — 35mm  and  y  =  35mm  and  are  assumed 
to  be  far  enough  away  from  the  sample  that  the  field  is  approximately  zero  at  these  boundaries. 
Indeed,  the  magnetic  flux  density  at  a  point  is  inversely  proportional  to  the  distance  between  the 
source  current  and  that  point.  Hence  the  field  tends  to  zero  as  the  distance  from  the  source  current 
or  conducting  sheet  increases.  A  schematic  of  the  resulting  two-dimensional  problem  is  depicted  in 
Figure  2  where  it  is  assumed  that  the  sample  (which  is  20mm  thick)  is  composed  of  aluminum  and 
the  conducting  sheet  (which  is  0.1mm  thick)  is  made  of  copper.  Thus  the  computational  domain 
which  we  will  use  for  the  purposes  of  developing  the  model  can  be  explicitly  defined  by 

=  {(x,  y,  z)  G  R3  :  0mm  <  x  <  50mm,  — 35mm  <  y  <  35mm}. 


-35mm 


Figure  2.  2-D  Schematic  of  Problem 

We  begin  the  development  of  the  model  by  introducing  a  mathematical  tool,  called  phasors, 
which  is  typically  used  in  the  held  of  electromagnetic  nondestructive  evaluation  whenever  periodic 
interrogating  inputs  are  employed.  As  mentioned  in  the  previous  section,  a  conducting  sheet  (copper 
in  our  example)  carrying  a  uniform  current  is  placed  above  the  sample  to  induce  eddy  currents  within 
the  sample.  Without  loss  of  generality,  we  assume  the  source  current  has  the  form 

Js  =  Jscos(uit)  k  =  JsRe(eluJt)  k 

where  Js  is  the  magnitude  of  the  source  current.  This  current  produces  a  magnetic  held  H(x,  y,  t ) 
described  by  Maxwell’s  equations.  As  the  magnetic  held  penetrates  into  the  sample,  a  phase  lag 
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occurs  due  to  the  finite  conductivity  of  the  sample  (aluminum  in  our  example).  In  other  words,  the 
magnetic  field  takes  the  form 

H(aj,  y,  t )  =  H(aj,  y)cos(ujt  +  0(x,  y)),  (1) 

where  H(x,  y)  is  a  vector  field  quantity  which  keeps  track  of  the  magnitude  and  direction  of  H  at 
each  point  in  space  while  9(x,y)  denotes  the  phase  shift  from  the  original  cosine  wave  at  the  same 
point  in  space  (this  term  takes  into  account  the  depth  of  penetration).  Since  these  are  the  quantities 
of  interest,  we  may  use  a  vector  phasor  ([1,  10]),  or  a  vector  of  complex  numbers,  H  defined  by 

fi(a;,  y,  t )  =  Re(H(x,  y)eiojt)  =  Im(U(x,  y)e^t+n^)  (2) 

which  keeps  track  of  only  these  quantities. 

We  then  use  Maxwell’s  equations  (they  are  the  basis  for  all  electromagnetic  phenomena)  in 
phasor  form  as  the  basis  for  our  derivation.  Details  on  the  derivation  of  Maxwell  equations  in 
phasor  form  from  Maxwell’s  equations  as  derived  from  first  principles  (e.g.,  Coulomb’s  law,  the 
Lorentz  transformation  and  relativity  theory  -  see  [14]),  can  be  found  in  [23].  We  then  have  for  the 
basis  of  all  the  following  derivations: 

V  •  B  =  0,  (3) 


V  •  D  =  p, 

V  x  E  =  —iuj B, 


(4) 

(5) 


and 


V  x  H  =  J  +  iuG.  (6) 

We  first  make  a  couple  of  remarks  regarding  (3)  -  (6).  To  begin  with,  our  system  is  considered  to 
be  electrically  neutral,  i.e.,  the  internal  electric  charge  density  p  equals  zero.  Secondly,  by  examining 
the  conductivity,  a,  of  aluminum  and  copper  (aai  =  3.72  x  107S'/m  and  ocu  =  5.8  x  10 7 S/m 
respectively)  and  by  using  Ohm’s  law 

J  =  ctE,  (7) 

we  can  argue  J  «  10 'E.  On  the  other  hand,  the  constitutive  law 

D  =  eE,  (8) 

where  e  is  the  electric  permittivity  (e  «  e0  «  ^  x  10_7F/m),  indicates  D  «  10_loE.  Therefore,  for 
the  source  frequencies  we  consider  within  the  scope  of  this  paper  (/<,  =  60Hz  -  2kHz),  uD  <  10_6E 
where  u  —  2nfs  is  the  angular  frequency.  Consequently,  in  the  sample  and  conducting  sheet 
J  >>  uD  which  implies  we  could  assume  uD  «  0  in  both  the  sample  and  conducting  sheet  in 
(6).  In  other  words,  the  term  uD  is  only  significant  in  the  air.  Both  the  literature  and  the  finite 
element  software  employed  in  the  computations  (Ansoft  Maxwell  2D  Field  Simulator)  neglect  the 
displacement  current  density  in  the  air  as  well  [28,  29,  30,  31,  47].  In  our  initial  computational  efforts, 
we  formulated  the  problem  both  including  as  well  as  neglecting  the  displacement  current  density 
and  compared  the  corresponding  solutions.  Our  findings  agreed  with  the  literature;  there  was  no 
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discernible  change  in  the  solution  when  the  displacement  current  density  was  ignored.  However,  the 
presence  of  this  term  does  alter  the  theoretical  results  slightly  (see  Section  3)  and  for  this  reason, 
we  choose  to  include  this  term  in  the  derivation  of  the  model. 

We  use  the  Ansoft  finite  element  solver  ([1])  in  our  computational  efforts;  therefore,  we  continue 
our  derivation  in  the  same  manner  as  done  in  Ansoft  by  introducing  a  magnetic  vector  potential  A. 
Based  upon  (3)  and  vector  null  identities,  B  can  be  represented  as  the  curl  of  a  vector  potential  A 
(called  the  magnetic  vector  potential), 

B  =  V  x  A.  (9) 

Therefore,  given  the  magnetic  vector  potential  A,  both  the  magnetic  held  H  and  magnetic  flux 
density  B  can  be  computed.  Accordingly,  we  want  to  combine  Maxwell’s  equations  to  obtain 
equations  in  conjunction  with  boundary  conditions  which  completely  determine  the  behavior  of  the 
magnetic  vector  potential  A  in  Q. 

Using  the  identity  B  =  V  x  A  in  (5),  we  have 

V  x  E  =  —  iu>(V  x  A)  or  V  x  (E  +  icoA)  —  0. 


Again,  using  vector  null  identities,  V  x  (E  +  iuiA)  =  0  implies  E  +  iuA  can  be  written  as  the 
gradient  of  a  scalar  potential,  denoted  by  4>.  As  a  result, 


E  =  —iuiA  —  V(f). 


(10) 


Finally,  we  can  use  (6)  and  (10)  in  conjunction  with  Ohm’s  law  (7),  the  constitutive  law  given 
by  (8)  and  the  constitutive  law  H  =  ^B  ( y  is  the  magnetic  permeability  in  H/m),  to  obtain 

V  x  V  x  A^j  =  (a  +  iu>e)(—iojA  —  \7(j))  Vx, y  e  0.  (11) 

In  the  above  equality,  the  right  side  represents  the  total  current  density  made  up  of  the  source 
current  density,  eddy  current  density,  and  displacement  current  density.  The  source  current  density 
Js  is  due  to  differences  in  electric  potential;  therefore,  Js  is  represented  by  the  term  —aWcj).  The 
term  —iua A  represents  the  eddy  current  density  Je  produced  by  a  time-varying  magnetic  held. 
Finally,  the  displacement  current  density  J a  due  to  time-varying  electric  helds  is  given  by  the  term 
iue(—iuA  —  Vd). 

Since  (11)  contains  two  unknowns,  A  and  (j>,  we  need  an  additional  equation  to  uniquely 
determine  solutions  of  the  system.  In  some  instances,  one  is  allowed  to  choose  a  “gauge”  which 
often  decouples  the  above  equation;  however,  based  upon  the  geometry  in  our  test  problem,  V- A  =  0 
is  naturally  imposed.  This  follows  since  the  only  nonzero  component  of  A  is  A3,  the  component 
of  A  in  the  z  direction  (the  direction  of  the  current  density  J).  Therefore,  V  •  A  =  —  0  by 

uniformity  in  the  z  direction.  Indeed,  this  is  the  Coulomb  gauge  [21,  pp.  221-222], 

Since  the  Coulomb  gauge  is  naturally  imposed,  it  provides  no  additional  information  in  our 
problem.  It  is  therefore  necessary  to  use  an  integral  constraint  to  obtain  a  second  equation  relating 
the  magnetic  vector  potential  A  to  the  scalar  potential  </>,  or  more  precisely  V</>.  For  the  integral 
constraint,  we  take  the  relationship 


I  —  J f  n da  —  /  (a  +  iue)(—iu A  —  V<? 6)  •  n da 

Js  Js 


(12) 
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where  /  is  the  total  measurable  current  flowing  in  a  conducting  surface  (S)  and  Jt  is  the  total 
current  density  within  the  conducting  surface.  Since  we  are  only  able  to  control  the  current  flowing 
in  the  conducting  sheet,  we  use  the  equivalent  equation  in  the  conducting  sheet 


Ts 


/  J<  •  n da—  /  (a  +  iuA  —  V</>)  •  nda 

J cs  J cs 


(13) 


where  Ics  is  the  total  current  flowing  in  the  conducting  sheet  (cs).  However,  since  (13)  only  provides 
a  relationship  between  A  and  in  the  conducting  sheet,  we  still  need  either  a  relationship  between 
the  potentials  or  a  condition  on  one  of  the  potentials  in  both  the  sample  and  air. 

Using  a  few  underlying  assumptions,  we  can  justify  V0  =  0  in  the  sample.  We  explore  two 
different  arguments  which  in  the  end  result  in  this  conclusion.  One  approach  is  to  examine  the  source 
current  density  term  caused  by  changes  in  potential  across  the  material  given  by  Js  =  —  ctV</>.  Since 
we  only  apply  a  current,  and  hence  effectively  a  potential,  across  the  conducting  sheet,  the  source 
current  density  (due  to  changes  in  potential)  only  exists  in  the  conducting  sheet.  Consequently, 
Js  =  —  ctV</>  =  0  in  the  sample.  Since  the  conductivity  of  aluminum  aai  does  not  equal  zero,  V(j) 
must  be  identically  zero  for  all  points  (x,  y )  in  the  sample.  A  second  alternative  is  to  assume  the 
sample  is  a  passive  conductor  modeling  a  short  circuit  (which  is  done  in  Ansoft).  In  general,  a 
passive  conductor  is  a  conductor  which  has  no  component  of  source  current  (as  discussed  above). 
In  other  words,  the  only  currents  considered  to  be  flowing  in  a  passive  conductor  are  eddy  currents 
and  displacement  currents  (when  considered).  In  addition,  a  short  circuit  conductor  is  treated  as 
a  conducting  ring  which  loops  back  on  itself.  In  this  case,  using  concepts  of  a  closed  circuit,  the 
change  in  potential  across  the  loop  is  zero  [1,  46],  i.e., 

V(f)  —  0  V(x,y)  G  sample.  (14) 

The  only  other  region  left  to  consider  is  the  air.  When  the  displacement  current  density  is 
ignored,  it  is  not  necessary  to  require  an  additional  constraint  on  V</)  in  the  air  since  (11)  reduces 
to 


V  x  | -V  x  A 


—  0  y(x,  y)  G  air. 


(15) 


However,  we  will  make  an  intuitive  argument  for  why  V<^>  should  be  taken  to  be  zero  in  the  air  as 
well.  We  use  an  argument  similar  to  that  one  used  for  the  sample.  The  source  current  density,  as 
discussed  previously,  is  only  present  in  the  conducting  sheet  and  hence  not  in  the  air.  However, 
by  examining  the  equations  alone,  this  does  not  give  us  any  knew  information.  It  simply  states 
that  Js  =  —aVcj)  must  equal  zero  in  the  air  which  is  already  true  since  a  is  identically  zero  in 
this  region.  On  the  other  hand,  we  can  intuitively  argue  that  the  source  current  density  is  due  to 
changes  in  potential,  and,  therefore,  if  no  source  current  density  is  present  in  a  region,  there  should 
be  no  change  in  potential  in  that  region  as  well.  Using  this  reasoning,  we  take  W(j)  —  0  in  the  air. 
Although  not  explicitly  stated  in  the  technical  notes  for  Ansoft  Maxwell  2D  Field  Simulator  ([1]),  the 
software  effectively  uses  V(j)  —  0  in  the  air  when  estimating  the  displacement  current  density  given 
the  approximate  finite  element  solution  A  determined  by  (13)  and  (14).  (In  the  actual  calculations 
produced  by  the  Ansoft  software,  displacement  current  density  is  ignored  totally.) 
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Therefore,  combining  all  the  above  we  have  two  coupled  equations  (11)  and  (13)  and  the 
additional  condition  V</>  =  0  for  points  ( x ,  y )  in  the  sample  and  air.  Using  these  equations,  the 
magnetic  vector  potential  A  can  be  uniquely  determined  if  appropriate  boundary  conditions  on  A 
are  specified. 

Recall  that  we  assume  evenly  symmetric  x  boundaries  due  to  the  symmetry  of  the  crack  and 
the  infinite  extent  of  the  materials.  In  other  words  on  the  x  boundaries,  we  assume  the  fields  on 
both  sides  of  the  boundary  oscillate  in  the  same  direction.  To  account  for  the  even  symmetry,  we 
assign  Neumann  boundary  conditions  to  these  boundaries.  In  a  similar  manner,  we  assume  the  y 
boundaries  are  “sufficiently  far”  away  from  the  sample  and  scanning  area  to  not  effect  the  overall 
measurements.  We  mentioned  previously  that  as  one  moves  farther  away  from  the  sample  and 
conducting  sheet,  the  magnetic  vector  potential  A  tends  to  zero.  Thus,  on  the  y  boundaries,  we 
assign  Dirichlet  boundary  conditions  to  indicate  the  boundary  is  “sufficiently  far”  away  from  the 
materials  so  that  A  fa  0.  Therefore,  the  magnetic  vector  potential  A  is  determined  according  to 

V  x  ^  J  V  x  A  (a:,  j/)^  =  {a(x,  y)  +  iue(x,  y))(-iuA(x,  y)  -  V</>)  x,y  G  U,  (16) 


J/  •  n  da 


/  (a(x,  y)  +  iue(x,  y))(-iuA(x,  y)  -  S7<fi) 

J  CS 


n  da 


and 


(17) 


with 


=  0  x,y  G  f2  \  cs 


A(x,  —35)  =  0 
VA-n|(0l„)  =  0 


=  A(x,  35) 

=  VA  •  n|(5o,„). 


(18) 


3.  Well-Posedness 


In  the  previous  section  we  developed  a  model  for  the  magnetic  vector  potential  A  given  a  source 
current  Ics.  The  boundary  value  problem  is  given  in  terms  of  two  unknowns  A  and  (f>  satisfying  (16) 
and  (17)  with  the  additional  condition  given  by  (18)  which,  as  we  shall  see,  uniquely  determines  the 
solution  when  coupled  with  appropriate  boundary  conditions.  However,  we  can  reduce  the  above 
system  into  a  single  integro-differential  equation  by  noting  (see  [23]  for  details)  that  the  term  V^>  is 
piecewise  constant  across  all  regions.  Therefore,  V</>  can  be  written  in  terms  of  A  in  (17).  Hence  we 
can  combine  all  the  equations  as  done  in  [23]  to  obtain  the  equivalent  form  of  the  boundary  value 
problem 

9  ~dx? ^  +  9  V ^ )  +  R(x,y)(a(x,y)  +  iuje(x,y))(iuMx,y)  +  K(x,y))  =  o,  (19) 

where  =  K(x,y)  is  defined  by 


K(x,y) 


_ Ics _ 

Acs  (acu+iue  cu) 


0 


for 

for 


{x,y)  e  cs 
(x,  y)  G  \  cs 


(20) 


Real  Time  Comp.  Alg.  for  Damage  Detection 


10 


Given  this  form  of  the  equation,  together  with  the  boundary  conditions,  we  consider  the 
existence  and  uniqueness  of  a  weak  solution  on  a  general  domain  given  by 

^  —  {("G  Vi  •  %min  Ci  %  Ci  %maxt  l/min  C  V  C  Umax} 

for  which  our  test  problem  is  a  specific  example.  Then,  let  H  =  L2(D)  and  V  =  {ip  G 
H1(Q)\'ip(x1  Umin)  =  0  =  ip(x,ymax)}  where  we  use  the  standard  Sobolev  space  notation,  if1(Q)  = 
{ip  G  L2(p)  :  Vip  G  L2(Q)}.  Note  that  we  interpret  pointwise  evaluation  of  functions  (along 
the  boundary  and  elsewhere)  in  terms  of  a  trace  operator  ([15])  for  which  we  suppress  notation 
throughout  this  paper  .  We  denote  by  (cp,  ip)  =  J2  (pipda  the  standard  inner  product  in  H  and 
(cp,ip)v  =  V0  •  \7ipda  the  (./^-equivalent)  inner  product  in  V . 

Then  using  integration  by  parts  together  with  natural  boundary  conditions  and  imposed 
conditions  on  test  functions  ip  EV,  the  variational  form  of  (19)  is  given  by 


,dA  dip .  ,dA  dip . \  ..  ,  .  . 

fa)  +  fyj)j  +  (luJT{o  +  *we)>l, ip)  +  (y(a  +  iuje)K , if)  =  0 

or  more  precisely 

(V^4,  Vip)  +  {PiA,  ip)  +  /32  /  ylda  /  ipda  =  /  fipda. 

J  cs  J  cs  J  cs 

where  ft  =  +  to*),  ft  =  alM]  f  =  if ala. 


(21) 


(22) 


3.1.  Existence  and  Uniqueness 

We  consider  the  existence  and  uniqueness  of  the  solution  ^4  to  (22)  (as  well  as  to  its  equivalent 
formulation  when  the  displacement  current  density  is  neglected)  in  the  context  of  a  Gelfand  triple 
setting  V  ^  H  ~  H*  V* .  We  have  that  the  embedding  V  H  is  dense  and  continuous  with 

| ip\H  C  k\ip\v  for  all  ip  £  V.  (23) 

where  the  norm  in  V  will  be  denoted  by  |  ■  ]y  and  |  •  |  will  denote  the  norm  in  H  for  the  rest  of  this 
section. 

We  define  a  sesquilinear  form  £  :  V  x  V  — »  (D  by 

£((p,ip)  =  (V(p,Vip)  +  {Pi<p,ip)  +  /32  f  (pda  (  ipda.  (24) 

J  cs  J  CS 

Then  (22)  can  be  written  as 

£(A,  ip)  =  f  fipda.  (25) 

J  CS 

We  intend  to  prove  £  is  G-continuous  and  for  certain  frequencies  is  G-elliptic.  We  will  then 
invoke  the  Lax-Milgram  theorem  to  prove  the  existence  and  uniqueness  of  the  weak  solution  A  to 
(22). 

Lemma  3.1  £  is  V -continuous,  i.e.,  there  exists  some  constant  c\  such  that 

\£{(p,ip)\  C  ci\(p\v\ip\v-  (26) 
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Proof:  Using  the  triangle  inequality,  Cauchy-Schwartz  inequality,  properties  of  L1,  Holder’s 

inequality  and  (23),  we  have 


£(<M) 


= 

(V^,  Vip) 

+  {fii (j>,  ip)  +  A  fcs  4>da  fcs  ipda\ 

< 

\(fi,'fi)v\  + 

\W 

i  A  ip)  1  +  1 A  fcs  fida  fcs  ipda\ 

< 

\fi\v\ip\v  + 

1 A 

oo\(PM  +  A  J'rs  ■<?  da  frs  ]tp  d 

< 

\C\v\*P\v  + 

1 A 

oo \<P  A  +  A ||0  nqn)  "0  nqn) 

< 

\PWVP\v  + 

IA 

oolAIA  +  |i|2IAII0lMl 

< 

(1  +  k2 A 

xi  T 

/c2jl|2|A|)|^|v|V’|v 

= 

ci\(p\v\ip\v, 

where  c\  —  1  +  /c2|A|<x  +  &2|1|2|A 


Lemma  3.2  There  exists  R  —  jF(/q  e,  Q)  such  that  for  fs  <  R(/j.,e,Cl),  there  exists  a  constant 
c2  >  0  such  that  C  satisfi.es 

WA)\  >cM2v.  (27) 

Proof:  By  (24)  and  the  definition  of  A  and  A 

ReC(<t>,  <f>)  =  (Vcj),  Vfi)  -  cv2(p,e(j),  </>)  +  u2f£f*cu  fcs  <fida  fcs  pda 


Therefore, 


Re£((f),(f))  >  (V0,  V<A  —  Lu2(iJ,e(f),  <p) 

>  (V^,V^)-a;2|/re|oo|^|2. 

Denoting  a  —  max(\ymin\,  \ymax\),  we  have  \y\  <  a  for  ( x ,  y)  €  Cl.  Hence  for  <p  e  V, 


\4>\2  —  JqI  ■  (jxpda 

=  [Xmax  fymax  1  •  (fibdydx 

J  xmin  ^  Vmin  _ 

=  -  Jn  VlFda  -  /ft  vfy4>do.  +  i;~ 

<  \(yf,A)\  +  \(y<P,%)\ 

<  2a|^||V0| 

<  ^!^|2  +  2a7!V^!2 

=  %\4>?  +  2a1\(t>\2v. 

Thus, 

(i  -  ^  \<t>?  <  2«7l0tv 


l^lv  _  _  l^l2  -  _  2a7)IAv- 


which  gives  us 
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1  -  7  = AM 


or 


2a 


7  = 


4(1  -  uo2 Ijueloo)' 
Note  that  the  condition  7  >  0  is  met  if  co  is  chosen  so  that 

1 

co  < 


In  order  to  satisfy  the  condition  7  <  A  an  initial  source  frequency  must  be  chosen  such  that 


co  < 


A -a2 


H 


or 


J_  1  -  a2 

S  <  27T  y  |/*e|oo 

where  co  =  2i r/s  and  we  tacitly  assume  a2  <  1.  Let  T  =  T (//,  e,  O)  =  then  for  /s  <  B , 

ReC(q J,  4>)  >  \(f\y  -  uo‘2\iie\oc\(f)\2  >  c2\(j)\y 


with  c2  =  1  —  2ay  >  0. 

We  note  that  for  our  test  problem,  p,  «  /i0  =  47r  x  10  t,  e  «  e0  ~  3^7  x  10-7  and  the  term  a 
discussed  in  the  proof  is  given  by  a  =  3.5  x  10~2.  Hence,  T  «  4.77  MHz.  All  of  our  computational 
results  employ  source  frequencies  fs  <  2 kHz.  Furthermore,  frequencies  much  larger  than  these  are 
not  appropriate  for  this  type  of  eddy  current  problem.  We  make  one  final  remark  about  the  bound 
T  on  source  frequencies  fs.  We  note  that  C  is  of  the  form  (( B  +  ip)v*,v'i  therefore,  a  sharper 
bound  on  the  source  frequencies  might  be  calculated  by  considering  the  spectrum  of  B  and  allowing 
those  frequencies  up  to  the  point  at  which  k  becomes  an  eigenvalue  of  B.  Since  the  bound  T  we 
have  calculated  is  adequate  for  analysis  including  all  calculations  we  consider  in  this  paper,  we  do 
not  choose  to  pursue  this  idea  here. 

A  consequence  of  the  Lax-Milgram  theorem,  as  discussed  in  [53,  pp.  271-275],  gives  an  existence 
and  uniqueness  theorem  for  solutions: 

Theorem  3.1  Let  V  H  V*  be  a  Gelfand  triple.  Let  f  xf-y( D  satisfy: 

(Al)  \C((j),  ip) |  <  ci\<l>\v\ip\v  for  all  (j),ip  eV, 

(A2)  \C(<t>,4>)\>  c2\f)\y  for  all  <t>  eV. 
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VO  =  <P,i>eV  (28) 

such  that  L  :  V  — V*  is  a  linear  topological  isomorphism  between  the  spaces  V  and  V* ,  and  for  the 
norms  we  have 


L\  <  Ci, 


--i 


1 

<  — . 

c2 


(29) 

(30) 


Put  otherwise,  the  weak  equation 

L<t>  =  /,  feV\ 

interpreted  as 

C(4>,  ip)  =  (/,  ip)  for  all  xp  e  V, 

possesses  for  each  f  G  V*  a  unique  solution  <p  G  V ,  and  this  solution  depends  continuously  on  f . 

Theorem  3.1  and  Lemmas  3.1  and  3.2  readily  give 

Theorem  3.2  There  exists  T  =  ^(/i,  e,  O)  such  that  for  fs  <  e,  Q),  there  exists  a  unique  weak 
solution  A  to  (22). 

We  make  a  note  that  when  the  displacement  current  density  is  neglected,  we  obtain  similar 
results  to  those  above  about  without  any  restrictions  on  the  source  current  density.  In  other  words, 
given  the  variational  form  (found  when  ignoring  the  displacement  current  density) 


(VA,Vip)  +  0  iA,iP)  +  fc 


/  Ada  /  ipda ,  =  / 

J  cs  J  cs  J  cs 


fipda 


(31) 


with  coefficients  /3i  =  iuuo  and  /32  =  —  w^ul7c“  we  obtain  the  following  theorem. 
Theorem  3.3  There  exists  a  unique  weak  solution  A  to  (31). 


3.2.  Continuous  Dependence  on  Parameters 

In  this  section,  we  give  theoretical  results  showing  the  solution  A  to  (22)  depends  continuously  on 
the  parameters  representing  the  damage  which  will  also  give  us  a  basis  for  convergence  results  in 
the  next  section. 

We  begin  by  considering  a  general  representation  of  a  damage  and  define  a  parameter  to 
represent  this  damage.  Although  we  only  consider  rectangular  damages  centered  along  the  length 
of  the  sample  in  our  computational  efforts,  we  will  allow  for  any  four  sided  polygon  or  quadrilateral 
in  the  theory.  We  will  then  relate  our  specific  problem  to  the  more  general  theory  considered  here. 

Any  quadrilateral  can  be  represented  by  its  four  corners.  Therefore,  we  let  q  = 
[(xi,  yi),  (x2, 1/2),  (^3, 2/3),  {xi,  2/4)]  be  a  vector  in  R4  x  R2  or  equivalently  R8  which  represents  the 
quadrilateral.  We  denote  by  Qad  the  set  of  admissible  parameters  q  where  it  is  assumed  Qad  is  a 
compact  subset  of  R8 .  We  note  that  it  is  possible  to  pick  Qad  to  be  compact  since  the  damage  is 
restricted  to  being  within  a  sample  of  finite  length  and  thickness.  We  also  allow  for  the  case  in 
which  no  damage  is  present  and  represent  this  case  by  “collapsing”  the  quadrilateral  into  a  single 
point  (with  zero  area),  i.e,  we  define  q  as  a  vector  of  identical  points. 
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Then  for  any  two  damages  given  by  parameters  q  and  q  in  Qad,  let 

d(q,q)  =  ||q-  q||  =  [(®i  -  ^l)2  +  (til  ~  Vi  f  +  •••  +  (®4  -  )2  +  (yA  -  yi)2]1/2  (32) 

to  be  the  standard  Euclidean  norm  in  R8.  We  denote  by  5Cl  (Figure  3)  the  points  in  Q  which  are 
either  in  the  damage  represented  by  q  or  q  but  which  are  not  in  both.  In  other  words,  let  fiq 
represent  the  points  (x,y)  in  Q  within  the  damage  given  by  q  and  flq  the  points  ( x,y )  in  Q  within 
the  damage  given  by  q.  Then  dCl  —  flq  U  flq  —  flq  fl  flq.  The  area  of  5Cl  is  directly  related  to 
the  distance  between  parameters  q  and  q.  Indeed,  d(q,  q)  — >■  0  if  and  only  if  x,\  — >  xi,  y\  — >-  yi, 
...,  Xi  —f  Xi,  %/i  — )-  yi,  i.e.,  the  four  corners  of  the  two  damages  approach  one  another.  Therefore, 
d(q,  q)  — >  0  implies  5 Cl  — >-  0  and  hence  the  area  of  5 Cl  — >  0. 


damage  q 

' — — n 

1 

damage  q 

♦  i 

Sample 

Figure  3.  The  Area  Represented  by 

To  relate  this  terminology  to  the  computational  examples  we  present  in  the  rest  of  this  paper, 
the  only  damages  we  consider  computationally  are  those  rectangular  in  shape  and  centered  along 
the  length  of  the  axis.  We  only  take  into  account  half  of  the  damage  such  that  the  left  boundary 
of  the  portion  we  consider  is  located  at  x  —  0  in  the  xy  plane.  Therefore,  the  only  variation  in  the 
damage  is  in  its  length  l,  thickness  h  and  depth  d. 

We  let  q  be  a  parameter  representing  only  these  variations,  i.e.,  q  is  a  vector  in  some  compact 
subset  Qad  of  either  R1,  R2,  or  R3.  For  example,  if  we  wish  to  only  estimate  the  length  of  the 
damage,  we  assume  the  depth  and  thickness  are  fixed  and  hence  q  =  q  —  l  is  in  a  compact  subset 
of  R1.  However,  if  we  wish  to  estimate  parameters  such  as  length  and  thickness  while  keeping  the 
depth  of  the  damage  fixed,  q  =  (l,  h )  and  hence  Qad  is  a  subset  of  R2.  If  we  allow  all  parameters  to 
vary,  q  =  (l,  h,d)  (E  Qad  C  R3.  Therefore,  given  q  (along  with  any  fixed  parameters)  we  have  l,  h, 
and  d  for  which  we  can  define  q  =  [(0,  —d),  (l,  —d),  (l,  — d  —  h ),  (0,  —  d  —  h)\  to  be  the  vector  of  corner 
points  for  the  damage.  This  yields  a  formulation  of  the  damage  in  terms  of  q  that  is  equivalent  to 
the  formulation  in  terms  of  q  given  with  metric  d.  Thus,  if  the  solution  depends  continuously  on 
q,  the  solution  will  also  depend  continuously  on  the  parameters  we  estimate  given  by  q. 

Let  V(5f2)  =  Hl(dCl)  with  the  pseudo  norm  =  J5^l'V(j)'V(f)da.  Then  given  the 

terminology  above,  we  first  prove  the  sesquilinear  form  C  depends  continuously  on  q  and  use 
this  result  to  obtain  the  continuous  dependence  of  A  on  q.  We  note  that  whether  we  include 
or  neglect  the  displacement  current  density  results  in  either  frequency  dependent  or  independent 
results,  respectively.  We  only  present  the  remaining  theory  with  the  displacement  current  density 
included. 
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Lemma  3.3  There  exists  a  positive  constant  K  independent  of  the  parameters  q,  q  such  that  for 
all  q,  q  G  Qad,  the  sesquilinear  form  £(q)(-:.  •)  satisfi.es 

|£(q)(<M)  -  £(q)(<M)l  <K\fi\  y(<5b)  Wv  (33) 

for  all  fi,  fi>  G  V  where  \  ■  |y^  — >  0  when  d( q,  q)  — >  0. 

Proof:  By  definition  of  the  sesquilinear  form  £  in  (24) 


|£(q)(^V0  -  ^(q)(^G 


|  (Vfi,Vfi)  +  (Pi{q}<j>,ip)  + 


/  fida  /  ipda\ 
J  cs  J  cs 


/  fida  / 

J  cs  J  cs 


■fida 


or  simply 


|£(q)(0,'0)  -  £(4)(<^)l  =  K(A(q)  -  A (4 ))</v 

Using  the  definition  of  /fi.  the  Cauchy-Schwartz  inequality  and  (23),  we  have 

K(A(4)  -  A(q))0,^)l  <  |(A(q)  - /3i(q))^||^| 

<  |*w||(Mq)(a(q)  +*we(q))  -  M4)(<K4)  +  *we(q)))^|  |-0| 

=  w\Tai(vai  +iU€al)  ~  *W^0e0||^|£2(in)l^l 

<  k2iw\piai(oai  +  icoeai)  -  iupl0e0\\fi\v{5Cl)\‘fi\v 

where  5Cl  is  defined  above.  Therefore,  if  we  let  K  =  k'2u\fj.ai(aai  +  iuitai)  —  iupL^efi,  we  have 

|((Al(4)  -  A(4 ))fi,fp)\  <  y(Jd)  IV’lv- 


(34) 


Theorem  3.4  Assume  the  admissible  parameter  set  Qad  is  a  compact  subset  of  R8.  Then  there 
exists  T  —  e,  U)  such  that  for  source  frequencies  fs  <  e,  Cl),  q  — >  A(q)  is  continuous  from 
Qad  to  1'  . 

Proof:  Let  q"-tqG  Qad  and  let  A(qn),  A(q)  be  the  corresponding  solutions  of  (25).  That  is, 
£(qn)04(qn),^)  =  for  ip  £  V  (35) 

£(4X^(4),  VO  =(/,'0)  for  ^Gb  (36) 

Subtracting  (36)  from  (35),  we  have 

^qn)(A(qn),fi)-C(q)(A(q),fi)  =  0. 

Adding  and  subtracting  the  term  £(qn)(A(q),  fi>)  and  simplifying,  we  have 

£(tf)(A(q")  -  A( q),  VO  =  -[£(4")  -  £(q)](A(q),  # 


Let  fi>  =  A{ qn)  —  A(q).  Then 

£(4n)(A((f)  -  A(q),  A(qn)  -  A( q))  =  -[£(4")  -  £(q)](^(q),  A(q")  -  A(q)). 
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Thus,  using  (33)  we  have 

|£(q")(^(q")  -  A(q),A(qn)  -  T(q))|  =  |[£(q")  -  £(q)](^(q),  H(q")  -  .4(q))| 

<  K\A(q) |y(5nn)|A(qn)  -  A( q)jy 

Then  we  can  use  the  above  equality  and  (27)  where  we  note  that  c2  can  be  established  independent 
of  qn  by  following  the  proof  of  Lemma  3.2  in  which  |//e|oo  can  be  taken  independent  of  qn.  Thus, 
defining  K  =  —  independent  of  qn,  we  obtain 

C2 

\A(qn)  -  -4(q)|y  <  K\A(cf)\v{gnn) 
where  the  right  side  goes  to  0  as  qn  — >  q. 

3.3.  Convergence  Results 

One  possible  estimation  problem  of  interest  consists  of  minimizing  over  some  set  Qad  of  admissible 
parameters  the  least  squares  functional 

n  m 

j(<i)  =  ~J2J2\A(Xi’yT(i)-Aij\2  (37) 

1 1  j= i 

where  {A13}  are  assumed  to  be  some  sampled  data  available  at  points  ( xi:yj):i  —  1  ,  ...,n, 
j  —  1  In  practice,  the  computations  for  the  minimizing  problem  are  carried  out  using 

an  approximate  system.  Here  we  will  consider  Galerkin  type  approximations  in  context  of  the 
variational  formulation  (22).  Let  HN  be  a  sequence  of  finite  dimensional  subspaces  of  H.  We 
denote  by  PN  the  orthogonal  projection  of  H  onto  HN .  Then  the  parameter  estimation  problem 
can  be  formulated  by  seeking  a  q  G  Qad  which  minimizes 

n  m 

jN (q)  =  o  J2 \aN(x^  vr  q)  -  Aij\2-  (38) 

i=  1  j= 1 

The  parameter  estimation  problem  given  above  assumes  we  have  sampled  data  A13  along  a  grid 
of  points  ( Xi ,  i/j),  i  =  1, ...,  n.  j  =  1, ...,  m  although  it  is  not  physically  possible  to  obtain  such  data. 
Typically  potentials,  for  example  A.  are  only  used  for  computational  purposes.  They  allow  one  to 
overcome  some  of  the  difficulties  which  arise  computationally  when  using  field  variables,  the  biggest 
complication  being  the  number  of  equations  which  must  be  solved  when  using  field  variables  [30]. 
The  field  B  or  H,  however,  is  the  only  measurable  quantity. 

Therefore  we  also  explore  an  alternative  problem  which  we  use  in  much  of  our  computational 
and  experimental  investigations.  It  involves  cost  functionals  using  observations  of  the  magnetic  flux 
density  B  =  V  x  A  =  (|^,  —  |^,  0)  =  (Lb,  B2,  0).  In  this  case,  (37)  and  (38)  are  replaced  by 


n  m 


1  n  iil 

j(h)  =  -  ^i2 

i-l  j- 1 

(39) 

1  n  m 

jN (q)  =  2  vr  q)  -  Bk  I2 

(40) 

i-1  j- 1 
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respectively,  for  k  =  1  or  2. 

However,  in  actual  experimental  situations,  the  sensor  detects  magnetic  flux  density  across 
the  entire  surface  area  of  the  sensor  and  then  outputs  this  data  as  data  for  a  specific  point.  This 
necessarily  involves  averaging  of  the  data  over  the  area  in  some  fashion.  Therefore  we  can  assume 

31*^1  B(x,y,cf)da  (41) 

j’jji  J&ij 

where  q*  represents  exact  parameter  values  corresponding  to  the  damage  within  the  sample  from 
which  the  data  was  obtained  and  Clij  is  a  neighborhood  of  ( x^yj ).  Hence  we  instead  consider  cost 
functionals  of  the  form 


n  m 


1  lb  lib 

i=l  j= 1 

(42) 

1  n  m 

jNm  =  2  EE|C'i-B"(q)--sb2 

(43) 

i= 1  3= 1 

as  well  as  the  corresponding  cost  functionals  in  terms  of  A  where 

CtJBk( q)  =  -J—  f  Bk(x,  y:  q )da. 

|“y|  dClij 

Thus,  using  the  notation  and  results  developed  in  this  section,  we  obtain  the  following  result. 
Theorem  3.5  Suppose 

(Bl)  The  finite  dimensional  subspaces  HN  satisfy  HN  C  V. 

(B2)  For  each  ip  G  V.  \PNip  —  ip\v  — *  0  as  N  — >  oo. 

Let  qjV  be  arbitrary  in  Qad  such  that  qjV  — >  q  in  Qad-  Then  there  exists  P  —  such  that 

for  source  frequencies  fs  <  e,  Cl),  we  have  Hw(qlV)  ->  H(q)  E  V  as  N  — oo,  where  ^4(q)  is  the 
solution  to  (25). 

Proof:  We  have 

£(q^)(^(q^),^)  =  (f,ipN)  for  ifN  G  HN  (44) 

£(q)(4f(q),-0)  =(/,'0)  for  ipeV.  (45) 

We  also  can  write 

\AN(qN)  -  A(q)\v  <  jHw(qw)  -  P*4(q)|y  +  \PNA(q)  -  A( q)|y. 

Hence  from  (B2)  it  suffices  to  prove 

TlV(qiV)  —  PNA(q)\v  — ■>  0  for  qjV  — )■  q  G  Qad  as  N  -»  oo. 

Taking  ip  =  ipN  in  (45)  and  subtracting  from  (44),  we  obtain 

£(qw)(H^(qw),^)  -  £(q)(^4(q),  ipN)  =  0 
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for  ipN  G  HN .  Furthermore,  by  adding  and  subtracting  £(qN)(PN  A(q),ipN)  and  £(qN)(A(q),ipN) 
and  simplifying,  we  have 

£(q^)(^(qw)  -  PNA(<i),r)  +  £(<lN)(PNA(<i)  -  A{ q)^N) 

+[£(q.N)-cmm^N)  =  o. 

Choosing  xfN  =  AN  =  AN(qN)  —  PNA(q),  we  hnd  that 

£(q^)(A7V,  AN )  +  £(qN)(PNA(q)  -  H(q),  AN )  +  [£(q")  -  £(q)](A(q),  AN)  =  0. 
Therefore, 

|£(qw)(Aw,  A")|  =  |  -  £(q")(P^4(q)  -  A(q),  AN)  -  [£(q")  -  £(q)](H(q),  AN)\. 

However,  using  the  triangle  inequality,  (26),  (33),  and  the  Cauchy-Schwartz  inequality,  we  can 
bound  the  right  side  of  the  above  by 

|  -  £(q")(PM(q)  -  A{ q),  AN)  -  [£(q")  -  £(q)](H(q),  AN)\  < 

(ci\PN A(q)  -  A(q)\v  +  K\A(q)\v{mN))\AN \v 

where  c\  is  the  constant  in  Lemma  3.1  and  K  is  the  constant  in  Lemma  3.3.  We  can  conclude  using 
(27)  that 

|A*V  <  -|P"^(q)  -  A(q)|v-  +  L|^(q)|„(Jfi„)|^(q)|l/. 

C2  C2 

Thus,  given  any  qiV  ->qG  Qad1  it  follows  from  (B2)  that  AN  =  AN(qN)  —  PNA(q)  — >  0  as  N  — »  oo, 
giving  us  the  desired  results. 

We  note  that  the  results  above  guarantee  convergence  of  both  A^q)  to  T(q)  in  V  and 
BAr(q)  ->•  B  (q)  in  L2(Q)  and  hence  provides  a  complete  theory  for  the  inverse  problem  considered 
(see  the  general  framework  given  in  Chapter  5  of  [5]). 

4.  Computational  Method 

Our  ultimate  goal  is  to  determine  the  feasibility  of  using  a  portable  sensing  device  in  conjunction 
with  inverse  problem  techniques  to  characterize  the  geometry  of  a  hidden,  i.e.,  subsurface,  damage 
within  a  sample.  To  achieve  this  goal,  we  must  develop  fast  and  efficient  forward  computational 
methods  to  be  used  possibly  numerous  times  in  the  inverse  problem  formulated  in  the  next 
two  sections.  To  this  end,  we  examine  reduced  order  Karhunen-Loeve  or  Proper  Orthogonal 
Decomposition  (POD)  techniques. 

The  POD  technique  is  an  attractive  order  reduction  method,  because  basis  elements  are  formed 
in  an  “optimal”  way  which  span  a  data  set  consisting  of  either  numerical  simulations  or  experimental 
data.  Since  the  POD  basis  is  formed  so  that  each  basis  element  captures  important  aspects  of  the 
data  set,  only  a  small  number  of  POD  basis  elements  are  needed  in  general  to  describe  the  solution 
[36].  Consequently,  if  the  POD  method  is  successful,  implementation  should  result  in  a  decrease  of 
computational  time. 


Real  Time  Comp.  Alg.  for  Damage  Detection 


19 


We  summarize  the  use  of  the  POD  method  in  the  context  of  the  least  squares  inverse  problems 
described  in  detail  in  Sections  5  and  6  and  first  introduced  in  [3,  4]  for  these  problems.  For  further 
details  on  the  general  POD  method,  we  refer  the  reader  to  [2,  6,  7,  12,  24,  26,  27,  32,  33,  34,  35,  36] 
and  the  extensive  list  of  references  contained  therein. 

The  first  step  in  forming  the  POD  basis  is  to  collect  “snapshots”  or  solutions  across  time, 
space  or  a  varied  parameter.  In  our  case,  we  let  q  be  the  vector  parameter  characterizing  physical 
properties  of  the  damage;  i.e.,  as  discussed  in  Section  3.2,  q  determines  the  geometry  of  the  damage 
including  the  length,  thickness,  depth,  etc.  of  the  damage.  For  an  ensemble  of  damages  {qj}^1; 
we  obtain  corresponding  solutions,  {A(qJ-)}^1,  of  (16)  with  (17)  and  (18),  for  magnetic  vector 
potentials  which  we  call  our  “snapshots”.  Alternatively,  from  the  solution  set  {A(qJ-)}-=!1,  we  can 
obtain  the  magnetic  fluxes  {B(qJ)}^1  and  instead  use  these  as  our  “snapshots”  if  we  wish  to  treat 
magnetic  fluxes  as  our  basic  state  variable.  (In  [3]  we  compare  results  using  one  field  versus  the  other 
as  the  basic  state  variable.  The  conclusions  are  summarized  in  Section  5.)  For  our  explanation,  we 
will  consider  snapshots  on  A  =  (0,  0,  A:i)  and  hence  our  explanation  will  be  for  the  scalar  case.  For 
the  vector  case,  we  would  simply  proceed  componentwise  [2,  12,  36].  Without  loss  of  generality,  we 
will  denote  the  vector  A  by  its  scalar  nonzero  component  A,  i.e.,  the  A3  component  of  A. 

We  seek  basis  elements  of  the  form 
Ns 

*i  =  E  v<  (46) 

3= 1 

where  the  coefficients  Vj(j)  are  chosen  such  that  each  POD  basis  element  Tj,  i  =  1,2,  ...,1VS, 
maximizes 

1  Ns 

K^qj);  ^)l2(o,  e)  I2 

S  3=1 

subject  to  ($j,  $j)i2(n,c)  =  1 1^*1 12  =  1-  Standard  arguments  guarantee  that  the  coefficients  Vj(j) 
can  be  found  by  solving  the  eigenvalue  problem 

CV  =  XV 

where  the  “covariance”  matrix  C  is  defined  by 

[C]  ij  =  -^-(A(qj),  A(qj))L2(n,c)- 

The  matrix  C  is  a  Hermitian  positive  semi-definite  matrix,  and  thus  it  possesses  a  complete 
set  of  orthogonal  eigenvectors  with  corresponding  nonnegative  real  eigenvalues.  We  order  the 
eigenvalues  along  with  their  corresponding  eigenvectors  such  that  the  eigenvalues  are  in  decreasing 
order, 

Ai  >  A2  >  ...  >  Ajvs  >0. 

We  furthermore  normalize  the  eigenvectors  corresponding  to  the  rule 
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Then  the  ith  POD  basis  element  is  defined  by  (46)  where  V](j)  represents  the  jth  component  of 
the  ith  eigenvector  of  C.  It  can  also  be  shown  that  {<hj}^1  are  orthonormal  in  L2(fi,(D)  and 
span^i}^  =  span{A(qj)}^v  Indeed,  given  any  A(qj)}  we  have 

Ns 

k= 1 


where 

ak(<lj)  =  (-4(qj),  Ta:)l2(o,c)- 

We  point  out  that  if  any  of  the  A*’s  are  zero,  say  A,  =  0  for  i  =  K  +  1, ...,  Ns,  then  even  though  the 
corresponding  V)  are  orthogonal  (and  of  course  linearly  independent),  we  will  have  span{$i}f=l  = 
span^i}^.  In  this  case  we  will  only  generate  K  <  Ns  linearly  independent  POD  basis  elements. 
A  discussion  of  the  relation  between  POD  basis  element  formation  and  the  popular  singular  value 
decomposition  (SVD)  methods  in  linear  algebraic  methods  is  given  in  [32], 

To  determine  the  reduced  number  N  of  POD  basis  elements  required  to  accurately  portray  the 
ensemble  of  “snapshots”  {A(qJ)}^1,  we  consider 

N  Ns 

£vl>  («) 

3= 1  j=1 

which  represents  the  percentage  of  “energy”  in  ■span{A(qJ)}^1  that  is  captured  in  span{Qj}f=1. 
The  reduced  basis  consists  of  only  the  first  N  elements  Tj,  i  =  1, ...,  N,  where  N  is  chosen  according 
to  the  percentage  “energy”  desired.  We  intuitively  argue  that  the  “energy”  we  are  referring  to  is 
related  to  the  total  electrostatic  energy.  Simply  stated,  the  matrix  C  contains  terms  of  the  form 

/  AAda  =  /  \A\2da. 

Jn  Jn 


which  can  be  written  in  terms  of  the  electric  field  E  according  to  (10),  E  =  —iuA  —  V(f).  Therefore, 


2da  —  Ci 


2da. 


where  C\  —  Since  V 6  is  piecewise  constant  (proved  in  [23]),  the  terms  in  the  matrix  C  are  a 
perturbation  of  terms  associated  with  electrostatic  energy  given  by 


We  —  -eo 


E\2dV. 


Thus  we  conclude  that  when  we  snapshot  on  the  magnetic  vector  potential,  the  ratio  in  (47)  is  a 
measure  of  the  electrostatic  energy  stored  across  D  (see  [14]). 

Employing  only  the  first  N  POD  basis  elements,  we  obtain  the  approximation  AN(cij)  for  A(q7) 
such  that 

N 

A( qj)  ~  ^"(qj)  = 

k-l 
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To  approximate  A^q)  where  q  is  a  given  parameter  not  in  the  set  {qj}^,  we  must  extend  the 
approximation  formula  to  obtain 

N 

^(q)  =  (48) 

k=l 


Two  possible  ways  of  computing  ak( q)  are  by  using  a  POD/Galerkin  method  or  a 
POD/Interpolation  method.  The  traditional  way  to  find  afc(q)  is  to  use  a  POD/Galerkin  method; 
however,  there  are  many  advantages  in  choosing  a  POD/Interpolation  method.  We  next  examine 
both  methods  and  compare  the  advantages  and  disadvantages  of  each. 

The  first  approach  we  consider,  the  POD/Galerkin  approach  is  essentially  an  application  of 
Galerkin’s  method  to  the  integro-differential  equation  (19)  in  conjunction  with  the  reduced  order 
POD  method.  The  POD/Galerkin  method  uses  the  approximation  given  in  (48)  in  the  variational 
form  of  the  integro-differential  equation  where  the  test  functions  are  chosen  to  be  the  reduced  order 
POD  basis  elements  -  The  system  then  reduces  to  a  linear  system  which  we  can  solve  for 

the  coefficients  a^q),  k  =  1 , ....  N. 

In  our  computational  efforts  reported  on  here  we  follow  the  literature  ([28,  29,  30,  31,  47]) 
and  neglect  the  displacement  current  in  the  numerical  implementation.  Therefore,  we  will  use  the 
variational  form 


1&4  dip  1  dA  dip  \  iu<jcu 

ad +  W  ad )  +  {  ^  ~  ^7 


f  Ada  I  ipda  =  -A-  f 

J  cs  J  CS  ^ cs  J  cs 


ipda. 


Substituting  (48)  into  (49)  and  letting  ip  =  <$*,  /  =  1, ...,  N,  we  obtain  the  system 


K  +  iuiM 


A , 


l-bb J 


a 


A  „ 


(49) 


(50) 


for  a  =  [«!,  «2, ...,  olnY  where 


and 


d$k  d$t  d$kd$u  , 
+  — - —)da  , 


^lk  (/2  /i  ^  dx  dx  '  dy  dy 

[M]ik  =  /  o$k$ida, 

Jn 

[b]i  =  /  &ida. 

J  CS 


Recall  p,  ~  //o  and  therefore,  changes  in  the  parameter  vector  q  only  change  the  conductivity  a  and 
therefore  only  effects  the  matrix  M.  Consequently,  in  the  inverse  problem,  for  each  “new  guess”  of 
q,  only  the  matrix  M  must  be  recalculated  with  each  iterative  step.  This  reduces  the  time  required 
for  each  forward  estimation  and  hence  the  total  time  for  the  entire  inverse  problem. 

Another  approach  to  forming  the  POD  approximation  is  to  use  POD/Interpolation  to  calculate 
the  coefficients.  This  method  relies  entirely  on  the  values  of  the  coefficients  ak( q)  for  q  in  the  set 
(qj}Si-  Unlike  the  POD/Galerkin  method,  it  does  not  take  into  account  the  boundary  value 
problem  which  A  satisfies. 
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Various  interpolation  methods  may  be  chosen  to  evaluate  cq,(q)  such  as  linear  interpolation, 
cubic  spline  interpolation  or  nearest  neighbor  interpolation.  However,  in  the  simulations  presented 
here,  the  built-in  Matlab  interpolation  functions  interpl  or  interp2,  one-dimensional  interpolation 
and  two-dimensional  interpolation  respectively,  were  used.  The  linear  interpolation  method  was 
chosen  for  the  one-parameter  simulated  case  and  cubic  spline  interpolation  was  chosen  for  all  other 
estimation  problems.  We  chose  these  methods  because  initial  trials  suggested  this  would  be  the 
best  method  to  choose  in  order  to  achieve  the  most  accurate  results  in  the  inverse  problem.  For 
a  more  detailed  discussion  on  these  interpolation  techniques,  we  will  refer  the  reader  to  [42,  pp. 
348-353]  and  [43,  pp.  93-106], 

A  detailed  discussion,  including  graphical  illustrations,  of  the  accuracy  of  the  POD 
approximations  compared  to  Ansoft  finite  element  approximations  using  each  method  can  be  found 
in  [23]  when  approximating  one  parameter  (length  of  damage).  Here  we  summarize  those  results. 
When  applying  the  POD/Galerkin  method,  using  only  N  =  3  POD  basis  elements,  we  were  able  to 
fairly  accurately  approximate  the  finite  element  solution  which  uses  over  7000  finite  elements  in  its 
approximation.  However,  as  the  value  of  N  increases,  the  approximation  continually  worsens  as  the 
conditioning  of  the  linear  system  being  solved  deteriorates.  When  N  =  3  basis  elements  are  used, 
the  condition  number  is  approximately  80;  however,  the  condition  number  jumps  to  approximately 
360  when  using  N  =  4  elements.  Moreover,  the  condition  number  is  over  12,000  when  all  of  the  POD 
basis  elements  (N  =  21)  are  used  in  the  approximation.  The  relative  error  in  the  approximation 
also  indicates  that  the  best  approximation  is  found  when  N  =  3  and  continually  worsens  for  larger 
values  of  N ;  therefore,  the  optimal  value  of  N  to  use  in  this  optimization  problem  appears  to  be 
N  =  3. 

The  POD/Interpolation  method  does  a  considerably  better  job  at  approximating  the  finite 
element  solution  with  only  2-3  basis  elements.  With  N  =  2  basis  elements,  there  is  still  some  visible 
error  in  the  approximation,  but  the  approximation  and  finite  element  simulation  are  approximately 
the  same,  with  less  than  a  1%  relative  error. 

From  the  examples  described  above,  there  are  cases  in  which  the  POD/Interpolation  method 
clearly  produces  more  accurate  results.  We  cannot  make  this  generalization  in  all  cases;  nonetheless, 
one  distinct  advantage  of  the  POD/Interpolation  method  is  that  it  does  not  rely  on  the  equations 
describing  the  system.  This  can  be  very  useful  in  some  experimental  applications  in  which  data  is 
available  but  it  is  not  easy  to  precisely  model  the  physical  process  corresponding  to  the  data.  In  this 
case,  if  there  is  correlation  in  the  data,  the  POD  method  may  be  a  viable  approximation  method  in 
which  an  appropriate  option  for  determining  the  coefficients  would  be  a  POD/Interpolation  method. 

5.  Simulated  Results 

In  this  section  we  present  computational  results  for  the  least  squares  inverse  problem  based  on 
the  methodology  developed  in  Section  4.  We  are  concerned  with  identifying  the  geometry  of  a 
crack,  i.e.,  estimating  parameters  such  as  the  length,  thickness,  and  depth  of  a  crack  within  a 
sample.  To  determine  the  feasibility  of  this  task  and  to  illustrate  the  use  of  the  reduced  order 
methodology,  we  first  estimate  a  single  parameter  keeping  the  other  two  parameters  fixed.  We 
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estimate  length,  thickness  and  depth  separately,  and  then  carry  out  a  two-dimensional  parameter 
estimation,  estimating  both  length  and  depth  simultaneously. 

In  [3]  and  [4],  we  performed  several  trials  in  which  we  assumed  we  had  access  to  various  types  of 
data,  such  as  the  A  held  or  the  B  held  at  various  points  {x^yf)  in  12.  We  compared  and  contrasted 
the  accuracy  to  which  we  could  estimate  the  length  l  of  the  damage  based  on  whether  the  A  held  or 
B  held  was  used  and  whether  we  considered  the  held  along  a  single  line,  multiple  lines  or  within  the 
entire  region  (which  is  not  physically  possible  and  was  only  tested  for  hypothetical  comparisons). 
From  the  results  in  these  references,  we  concluded  that  extremely  accurate  results  were  obtained 
only  when  the  y  component  B2  of  the  magnetic  hux  density  was  used  in  the  cost  criterion,  i.e., 
when  we  used 

1  n  m 

J(q)  =  io8£f  (^,  %,  q)  -  io8B2J(q*)|2  (51) 

i= 1  3= 1 

where  108  is  a  scaling  factor  accounting  for  the  low  order  of  magnitude  of  the  held  (B-2  is  on  the 
order  of  10_8lF6/m),  B %  (q)  is  the  reduced  order  POD  approximation  to  the  magnetic  hux  density 
given  by  BN  —  V  x  AN ,  and  B  is  “data”  from  a  sample  we  wish  to  characterize.  (In  this  section, 
B  is  obtained  from  finite  element  simulations  with  q  =  q*  to  which  randomly  generated  noise  has 
been  added  in  the  usual  manner  (see  [4]  or  [23]).  Furthermore,  performing  multi-line  scans  or  using 
full  region  data  improved  the  results  only  marginally  and  hence  did  not  warrant  the  extra  effort 
and  time  in  collecting  more  extensive  data  sets.  Consequently  the  results  presented  in  this  section 
involve  only  the  least  squares  difference  in  the  B2  held  given  by  (51)  along  a  single  line  located 
1mm  above  the  conducting  sheet. 

We  discussed  above  two  different  methods  used  in  forming  the  reduced  order  POD 
approximation:  the  POD/Galerkin  method  and  POD/Interpolation  method.  As  mentioned 
previously  the  POD/Interpolation  method  did  much  better  at  approximating  the  finite  element 
solution;  therefore,  we  only  present  here  the  inverse  problem  results  obtained  using  that  method. 
We  refer  the  interested  reader  to  [23]  for  a  summary  of  and  comparison  with  results  using  the 
POD/Galerkin  method. 

5.1.  Estimating  One  Parameter 

In  the  one-parameter  estimation  problem,  we  estimated  three  different  lengths  (T  =  1.3mm, 
T  =  2.5mm  and  T  =  5mm),  one  thickness  (h*  =  1.3mm)  and  three  different  depths  (d*  =  3mm, 
d*  —  8mm  and  d*  —  11mm)  while  keeping  all  other  parameters  fixed.  In  estimating  the  length 
of  the  damage,  we  generated  an  ensemble  of  damages  keeping  the  thickness  fixed  at  2 mm  with 
various  crack  lengths  at  a  fixed  depth  of  9mm.  We  varied  the  lengths  from  Omm  to  4mm 

in  increments  of  0.2mm  resulting  in  Ns  —  21  different  damages  and  use  the  commercial  software 
Ansoft  Maxwell  2D  Field  Simulator  to  generate  snapshots  {A(/j)}^21.  Similarly,  in  estimating  the 
thickness  of  the  damage,  we  kept  the  length  fixed  at  2 mm  and  depth  fixed  at  9mm  and  varied  the 
thickness  from  Omm  to  4mm  in  increments  of  0.2mm.  Finally,  in  estimating  the  depth,  we  kept 
the  length  of  the  damage  fixed  at  1.5mm  and  thickness  fixed  at  0.5mm.  We  took  snapshots  on  A 
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for  damages  with  depths  ranging  from  0.25mm  to  19.25mm  in  increments  of  0.5mm  (Ns  =  39  total 
snapshots).  Since  we  added  random  error  to  the  data  to  mimic  normal  random  measurement  error 
found  in  experimental  data,  we  generated  ten  different  sets  of  random  error  which  was  added  to  the 
data  and  averaged  the  results  across  the  ten  trials  to  obtain  an  average  result  along  with  a  standard 
deviation  across  the  ten  trials.  (See  [23]  for  full  details.)  Using  the  optimization  Nelder-Mead 
([25,  41]),  we  obtained  the  results  in  Table  1. 


Table  1.  Results  Estimating  a  Single  Parameter  Keeping  All  Other  Parameters  Fixed  at  a  10% 
Relative  Noise  Level 


q 

(in  mm) 

No.  POD 

True  q* 
in  data 

Avg.  Est. 

q 

St.  Dev. 

length 

4 

1.3 

1.2977 

0.0057 

4 

2.5 

2.4981 

0.0020 

4 

5.0 

4.9771 

0.0066 

thickness 

9 

1.3 

1.3056 

0.0054 

depth 

5 

3.0 

3.0349 

0.0115 

5 

8.0 

8.0631 

0.0109 

5 

11.0 

10.9184 

0.0107 

We  note  that  in  each  case,  we  considered  parameters  not  included  in  those  upon  which  we 
snapshot.  Furthermore,  we  consider  one  case  in  which  the  parameter  value  was  outside  the 
interval  upon  which  we  took  snapshots.  In  other  words,  when  estimating  length  T  =  5mm, 
the  POD/Interpolation  method  could  not  be  used.  We  instead  considered  POD/Extrapolation. 
Nonetheless  in  all  the  results  presented,  the  method  was  shown  to  be  quite  accurate  with  low 
standard  deviations  in  all  cases.  We  also  note  that  although  we  have  reported  on  three  estimated 
depths  in  which  we  obtained  extremely  accurate  estimates,  in  some  of  the  cases  tested,  the  results 
were  not  quite  as  accurate.  For  all  the  depths  tested  ranging  from  d*  =  1mm  to  d*  =  8mm 
we  obtained  accurate  results.  For  depths  past  8mm,  there  was  no  readily  identifiable  criterion  to 
predict  a  priori  which  depths  could  be  estimated.  Some  estimates  of  d*  in  the  range  9mm  to  20mm 
were  good  and  some  estimates  would  be  considered  only  fair  to  poor.  Full  details  can  be  found  in 
[23].  Nonetheless,  the  results  indicate  the  viability  of  this  method  in  single  parameter  estimation 
problems. 

5.2.  Estimating  Two  Geometric  Parameters 

In  Section  6  we  discuss  a  need  to  modify  the  assumptions  made  in  the  original  test  problem  to 
more  accurately  describe  the  behavior  of  experimental  data  obtained.  In  short,  the  computational 
domain  was  expanded  beyond  the  edges  of  the  sample  and  snapshots  were  taken  of  the  magnetic 
flux  density  data  only  on  a  single  line  above  the  conducting  sheet  (instead  of  the  whole  region).  For 
the  computational  trials  in  the  1-D  case,  we  took  snapshots  of  the  magnetic  vector  potential  for  the 


Real  Time  Comp.  Alg.  for  Damage  Detection 


25 


entire  computational  domain  even  though  in  the  inverse  problems  we  only  used  those  data  points 
along  a  single  line.  Furthermore,  in  dealing  with  experimental  data  we  considered  data  across  the 
entire  length  of  the  sample,  instead  of  just  half  the  sample  as  done  in  1-D  computational  examples 
previously.  We  also  implemented  these  changes  in  the  two-parameter  estimation  problems. 

Otherwise,  we  proceed  as  in  the  previous  estimation  problems  by  first  generating  an  ensemble 
of  damages.  We  consider  damages  with  depths  ranging  from  1  mm  to  11mm  in  increments  of  2mm 
in  combination  with  lengths  from  0.5cm  to  3.5 cm  in  increments  of  1cm  (we  now  consider  longer 
damages  similar  to  those  in  the  next  section  involving  experimental  data).  We  keep  the  thickness 
fixed  at  1mm.  A  total  of  24  snapshots,  {B2(di,lj)},  i  =  1, ...,  6,  j  =  1, ...,  4  were  generated  using 
Ansoft.  We  present  results  in  Table  2  for  estimating  a  depth  of  d*  =  2 mm  and  length  T  =  1  cm,  a 
depth  of  d*  =  4mm  and  length  T  =  2 cm,  and  finally  a  depth  of  d*  =  6mm  and  length  T  =  3 cm. 
(In  each  case,  10%  relative  noise  was  added  to  the  generated  data  before  use  in  the  inverse  problem 
criteria.) 


Table  2.  Results  Estimating  Two  Parameters  Simultaneously  at  a  10%  Relative  Noise  Level 


q 

No.  POD 

True  q* 

Avg.  Est.  q 

St.  Dev. 

depth 

5 

2  mm 

2.0473mm 

0.0045mm 

length 

5 

1  cm 

1.0180cm 

0.0087cm 

depth 

5 

4  mm 

4.0850mm 

0.0074mm 

length 

5 

2  cm 

1.9801cm 

0.0098cm 

depth 

5 

6  mm 

6.0316  mm 

0.0031mm 

length 

5 

3  cm 

3.0396cm 

0.0055cm 

5. 3.  Conclusions 

When  using  the  B2  field  in  the  inverse  problem,  the  methods  proved  to  be  accurate  and  robust, 
allowing  us  to  accurately  estimate  the  length,  thickness,  and  depth  of  a  damage  within  a  sample  as 
well  as  length  and  depth  simultaneously  even  when  the  data  contained  considerable  noise. 

Furthermore,  there  are  two  significant  findings  to  report  with  regard  to  the  above  results. 
First  of  all,  in  most  cases  we  were  able  to  use  10  POD  basis  elements  or  less  in  each  of  the  trials 
performed.  We  compare  this  to  a  total  of  over  7000  finite  elements  required  to  solve  the  boundary 
value  problem  initially  (using  Ansoft  Maxwell  2D  Field  Simulator).  Hence,  if  one  were  to  use  the 
finite  element  software  for  each  forward  run,  we  could  expect  a  time-intensive  inverse  problem. 
This  leads  us  to  the  second,  and  most  signifi,canf  finding  with  regard  to  reduction  in  computational 
time  which  can  be  summarized  as  follows.  If  one  were  to  use  a  software  package  such  as  Ansoft ’s 
Maxwell  2D  Field  Simulator  to  calculate  the  forward  problem  each  time  it  is  required  in  the  inverse 
problem,  it  would  take  approximately  5-7  minutes  for  a  single  forward  solve  and  hence  any  inverse 
algorithm  based  on  this  forward  solver  may  require  several  minutes  to  several  hours  of  time  for  the 
optimization  problem.  In  using  the  reduced  order  POD  methodology  for  the  forward  problem,  the 
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entire  inverse  problem  takes  approximately  8  seconds,  less  than  T  the  time  required  for  a  single 
forward  simulation.  As  a  forward  algorithm  is  called  numerous  times,  this  is  a  substantial  reduction 
in  time  required.  For  example,  assuming  an  average  of  20  iterative  steps  in  the  typical  optimization 
procedure  for  these  problems  and  a  total  forward  simulation  time  of  5  minutes  (300  seconds)  to 
7  minutes  (420  seconds)  for  a  finite  element  forward  solve,  the  total  computational  time  for  the 
inverse  problem  would  range  from  1  hour  40  minutes  (6000  seconds)  to  2  hours  20  minutes  (8400 
seconds).  Thus,  we  arrive  at  a  speed  up  factor  ranging  from  750  to  1050,  a  factor  of  approximately 
103. 

Furthermore,  most  of  the  extensive  computational  time  is  required  only  in  the  initial  collection 
of  snapshots  which  would  take  place  prior  to  implementation  in  a  practical  setting.  This  suggests 
that  a  portable  sensing  device,  when  coupled  with  reduced  order  modeling  in  the  inverse  problem, 
might  be  plausible  in  practical  damage  detection  applications. 

6.  Experimental  Results 

The  simulations  performed  in  the  previous  section  are  most  encouraging  and  suggest  a  natural  next 
step  to  further  test  our  proposed  methods  with  experimental  data.  Therefore,  we  designed  the 
experiment  depicted  in  Figure  4,  in  which  we  try  to  detect  and  parameterize  a  damage  within  an 
aluminum  sample  using  a  giant  magnetoresistive  (GMR)  sensor.  The  sample  is  constructed  of  17 


Figure  4.  Experimental  Setup 

layers  of  1  mm  thick  aluminum  plates  with  a  slice  cut  out  of  one  of  the  layers  to  simulate  a  damage 
within  the  sample  (see  Figure  5).  The  “damaged”  piece  of  aluminum  is  moved  from  one  layer  to 
another  to  simulate  damages  within  the  sample  at  different  depths,  and  the  length  of  the  damage  is 
varied  by  producing  “gaps”  of  varying  size  from  the  aluminum  plate  (the  thickness  of  the  damage 
is  always  fixed  at  1mm).  As  a  means  of  inducing  current  within  the  sample,  a  thin  sheet  of  copper 
carrying  a  uniform  current  of  3 A  is  placed  above  the  sample  on  top  of  a  thin  sheet  of  paper  (to  avoid 
direct  physical  contact  between  the  sample  and  the  conducting  sheet).  The  GMR  sensor  measures 
the  amplitude  and  phase  of  the  magnetic  flux  density  across  a  2 in  line  (along  the  length  of  the 
sample)  every  0.635mm.  The  data  is  then  filtered  through  a  lock-in  amplifier  and  saved  to  a  file. 
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Furthermore,  we  took  data  across  various  frequencies  to  analyze  the  effect  of  the  frequency  on  the 
estimation  problem. 


Figure  5.  Schematic  of  the  Damaged  Layer 

Frequency  is  a  very  important  factor  to  consider  when  trying  to  detect  and  parameterize  a 
damage.  This  is  due  to  the  depth  of  penetration  of  the  eddy  currents,  also  referred  to  as  the 
“skin  effect” .  Eddy  currents  are  not  uniformly  distributed  throughout  a  material  but  instead  decay 
exponentially  with  depth  in  the  material  [40].  The  distance  at  which  the  eddy  current  density  has 
decreased  by  a  factor  of  1/e  (36.8%)  is  called  the  depth  of  penetration  and  can  be  calculated  by 

e  =  -fA=  (52) 

VnfsRC 

where  fs  is  the  source  frequency,  /i  is  the  magnetic  permeability  and  a  is  the  conductivity  of  the 
material  [10,  p.  370].  We  carried  out  experiments  using  four  different  frequencies,  250Hz,  500Hz, 
1kHz  and  2kHz  with  a  depth  of  penetration  ranging  from  2.70 mm  to  7.64mm  for  2kHz  down  to 
250Hz,  respectively.  This  provides  varying  sensitivity  when  considering  damages  at  different  depths. 

Given  the  magnetic  flux  density  data,  from  the  GMR  sensor  for  a  given  damage  at  a 
specified  depth  d*  and  with  a  given  length  T,  we  wished  to  estimate  these  parameters  using  an 
appropriate  cost  criterion.  The  first  step  in  the  optimization  process  was  to  generate  snapshots 
representative  of  the  experimental  data  across  the  various  damages.  To  generate  the  snapshots,  we 
first  explored  the  idea  of  using  simulations  obtained  from  the  finite  element  solver  Ansoft  Maxwell 
2D  Field  Simulator  to  form  the  POD  basis  elements  as  done  in  the  previous  section.  However,  in 
order  for  the  snapshots  to  be  representative  of  the  data,  we  needed  to  modify  the  assumptions  made 
in  the  original  test  problem.  We  had  originally  assumed  a  sample  of  infinite  length  to  disregard 
boundary  effects  from  the  edges  of  the  material.  Unfortunately,  the  experimental  data  showed 
significant  boundary  effects.  Therefore,  we  modified  the  computational  domain  used  in  the  finite 
element  solver  to  that  domain  depicted  in  Figure  6  which  includes  the  edges  of  the  sample.  Moreover, 
instead  of  only  considering  half  of  the  sample,  we  now  considered  data  across  the  entire  length  of 
the  sample  as  collected  in  the  experimental  setup. 

However,  after  analyzing  the  experimental  data  and  the  Ansoft  simulations,  we  noticed 
significant  differences  between  the  data  and  simulations.  (For  details  and  analysis  of  these 
differences,  see  [23].)  Consequently,  we  chose  to  instead  snapshot  directly  on  the  experimental 
data  to  form  the  POD  basis  elements.  Furthermore,  to  obtain  a  definite  pattern  in  the  data  as 
a  function  of  the  damage  within  the  sample,  it  was  necessary  to  filter  out  the  background  noise 
(data  obtained  when  the  sample  contained  no  damage)  and  to  average  the  amplitude  across  the 
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Figure  6.  Modified  Computational  Domain 


center  of  the  damage  (this  produces  desired  symmetry  in  the  data).  However,  while  filtering  out 
the  background  noise  enhanced  the  differences  in  the  data  as  a  function  of  the  various  damages,  at 
the  same  time,  it  also  seemed  to  intensify  the  boundary  effects  present.  Thus,  it  was  necessary  to 
disregard  the  data  at  the  edges  of  the  sample.  In  other  words,  we  consider  the  cost  function  given 
by 


j=a- 2 


(53) 


where  q  is  the  vector  containing  the  parameters  we  wish  to  estimate,  Bff  (q)  is  the  POD 
approximation  formed  using  snapshots  on  the  data  itself,  B \  is  GMR  data  at  grid  points  x3, 
j  =  1, ...,  n  with  n  total  grid  points  and  a  and  b  indicate  the  indices  of  the  grid  points  we  consider 
in  our  cost  criterion. 

In  [23],  we  also  analyzed  the  data  when  the  phase  data  was  manipulated  with  either  a  phase  shift 
and/or  when  the  phase  data  was  averaged.  We  took  data  at  different  frequencies  (the  importance 
of  which  was  mentioned  above),  and  depending  on  the  frequency,  either  the  real  portion  of  the  data 
or  imaginary  portion  of  the  data  exhibited  sporadic  behavior.  This  was  due,  in  part,  to  the  phase 
of  the  data.  If  an  appropriate  phase  shift  was  incorporated,  the  sporadic  behavior  of  the  data  could 
be  minimized.  Furthermore,  since  it  was  necessary  to  average  the  amplitude  across  the  center  of 
the  damage,  it  seemed  logical  to  also  consider  averaging  the  phase  data  across  the  center  of  the 
damage.  We  performed  extensive  analysis  of  the  data  with  these  various  modifications  in  [23]  and 
concluded  that  the  best  results  of  the  inverse  problem,  in  all  but  one  case  considered  (the  case  in 
which  a  source  frequency  of  1kHz  was  used),  could  be  found  when  the  data  incorporated  a  phase 
shift.  The  phase  shift  seemed  to  intensify  the  variations  in  the  data  as  a  function  of  the  damage 
within  the  sample.  Furthermore,  it  was  necessary  to  use  only  the  real  portion  of  the  data  or  the 
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imaginary  portion  of  the  data  in  the  cost  criterion,  but  not  both  parts  (due  to  the  sporadic  nature 
of  the  data  as  discussed  above).  The  choice  of  either  the  real  portion  or  imaginary  portion  could 
be  readily  determined  by  simply  examining  the  data  before  implementation  and  determining  which 
portion  exhibited  a  clear  pattern  as  opposed  to  sporadic  behavior  which  would  be  found  in  the 
other  portion  of  the  data.  In  the  case  of  1kHz,  a  phase  shift  did  not  intensify  the  variation  in  the 
data  as  in  the  other  cases  but  did  just  the  opposite.  As  a  consequence,  for  this  frequency  it  was 
necessary  to  consider  the  data  with  no  phase  modifications. 

Before  presenting  results  for  the  inverse  problems,  we  make  one  final  note  regarding  the  use  of 
only  a  portion  of  the  data  in  the  cost  criterion.  To  apply  the  results  in  a  practical  setting,  it  would 
be  necessary  to  either  develop  a  model  which  accurately  portrays  the  data  or  to  collect  data  over 
a  period  of  time,  building  a  data  base  which  inherently  contained  the  variations  in  the  data.  With 
this  extensive  data  set,  we  would  be  able  to  discern  the  pattern  in  the  data  a  priori.  Thus,  the 
frequency  would  indicate  the  appropriate  form  of  the  data  to  be  used  in  the  cost  criterion. 

Next,  we  present  representative  results  below  for  the  various  frequencies,  in  which  we  use  only 
the  real  portion  of  the  data  with  a  phase  shift  incorporated  in  the  cost  function  given  by  (53)  for 
frequencies  250Hz  and  500Hz,  the  imaginary  portion  of  the  data  with  no  phase  shift  is  used  in  the 
cost  function  for  frequency  1kHz  and  the  imaginary  portion  of  the  data  with  a  phase  shift  is  used 
in  the  cost  function  for  frequency  2kHz. 

6.1.  Determining  the  Length  of  the  Damage 

In  determining  the  length  of  the  damage,  we  kept  the  depth  fixed  at  2 mm  or  3 mm  and  estimated 
a  length  of  either  lcm  or  1.5cm  (keeping  the  thickness  fixed  at  1mm).  When  estimating  length, 
snapshots  for  the  POD  basis  elements  incorporated  data  from  samples  with  damages  at  the  fixed 
depth  with  varying  lengths,  excluding  the  true  length.  In  other  words,  in  trying  to  estimate  a 
length  of  lcm  at  a  fixed  depth  of  2 mm,  we  used  data  from  samples  with  damages  having  lengths  of 
0.5cm,  1.5cm,  2 cm  and  3 cm  all  at  a  depth  of  2 mm  to  form  the  snapshots.  In  addition,  all  the  POD 
basis  elements  were  used  in  the  approximation  since  we  had  only  4  basis  elements.  (Only  a  small 
amount  of  data  was  taken  to  establish  proof-of-concept.)  The  results  of  the  estimation  problem 
can  be  found  in  Table  3  (where  indicates  that  the  optimization  routine  failed  to  converge  to  an 
estimated  parameter). 

6.2.  Determining  the  Depth  of  the  Damage 

In  determining  the  depth  of  the  damage,  we  analyzed  the  results  obtained  when  detecting  the  same 
depth  using  various  frequencies  at  a  fixed  length  of  lcm  and  1.5cm.  In  forming  the  snapshots, 
we  did  the  same  as  with  length  where  our  snapshots  were  given  by  magnetic  flux  density  data  for 
samples  with  damages  with  the  chosen  fixed  length  (and  fixed  thickness  of  1mm)  at  all  the  depths 
(1mm,  2 mm,  3mm,  4mm,  6mm,  and  8mm)  excluding  the  depth  we  wish  to  estimate.  Tables  4  and 
5  summarize  the  results. 
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Table  3.  Determination  of  Length 


Frequency 
(Hz)  ' 

Depth 

(mm) 

Actual  Length 
(cm) 

Optimized  Length 

(cm) 

Relative  Error 

250 

2 

1.0 

1.5068 

50.68% 

1.5 

2.5713 

71.42% 

3 

1.0 

0.9254 

7.46% 

1.5 

1.3869 

7.54% 

500 

2 

1.0 

0.8002 

19.98% 

1.5 

1.6036 

6.91% 

3 

1.0 

0.9225 

7.75% 

1.5 

1.5540 

3.60% 

1000 

2 

1.0 

0.8169 

18.31% 

1.5 

1.4789 

1.41% 

3 

1.0 

- 

- 

1.5 

- 

- 

2000 

2 

1.0 

0.7566 

24.34% 

1.5 

2.7050 

80.34  % 

3 

1.0 

1.3782 

37.82% 

1.5 

1.3100 

12.67  % 

Table  4.  Determination  of  Depth  with  Fixed  Length  1.0cm 


Actual  Depth  (mm) 

Frequency  (Hz) 

Optimized  Depth  (mm) 

Relative  Error 

2 

250 

0.9411 

52.95% 

500 

2.1919 

9.59% 

1000 

2.1191 

5.96  % 

2000 

2.0479 

2.39  % 

3 

250 

3.4827 

16.09% 

500 

2.8047 

6.51  % 

1000 

2.9004 

3.32  % 

2000 

2.9127 

0.91  % 

6. 3.  Determining  the  Length  and  Depth  of  the  Damage 

In  our  final  trials,  we  estimated  both  length  and  depth  simultaneously.  In  order  to  implement  the 
interpolation  routine  for  our  POD  approximation  of  the  magnetic  flux  density,  we  need  to  have  a 
grid  upon  which  we  know  the  values  of  our  coefficients  for  our  POD  approximations.  However,  in 
forming  our  snapshots  for  a  specific  estimation,  we  exclude  all  the  data  with  either  the  same  depth 
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Table  5.  Determination  of  Depth  with  Fixed  Length  1.5cm 


Actual  Depth  (mm) 

Frequency  (Hz) 

Optimized  Depth  (mm) 

Relative  Error 

2 

250 

1.5721 

21.40% 

500 

1.8979 

5.10  % 

1000 

1.8330 

9.35  % 

2000 

1.6484 

17.58% 

3 

250 

- 

- 

500 

3.1094 

3.15% 

1000 

3.1274 

4.25% 

2000 

3.2070 

6.90% 

or  length  as  that  we  are  attempting  to  estimate.  For  example,  if  we  are  estimating  a  length  of 
T  —  1.5 cm  and  a  depth  of  d*  —  2 mm,  our  snapshots  consist  of  data  for  samples  with  damages 
having  lengths  0.5cm,  1cm,  2 cm  and  3 cm  at  depths  1mm,  3mm,  4mm,  6mm  and  8mm  (all  with 
fixed  thickness  of  1mm).  In  this  case,  we  only  analyze  the  results  using  a  frequency  of  500Hz.  The 
results  can  be  found  in  Table  6. 

Table  6.  Determination  of  Depth  and  Length  Simultaneously  with  Frequency  500Hz 


Actual  Length  l 

1cm 

1.5cm 

2mm 

l  =  1.0635 

l  =  1.8080 

Actual 

d  =  2.3097 

d  =  1.8403 

Depth  d 

3mm 

l  =  0.9065 

d  =  2.9522 

l  =  1.4612 

d  =  2.9759 

6-4-  Conclusions 

Assuming  a  fixed  source  frequency  and  its  associated  cost  criterion,  we  were  able  to  demonstrate 
that  the  POD  method  in  the  context  of  inverse  problems  is  a  viable  method  even  with  experimental 
data.  Depending  on  the  frequency  used,  we  were  able  to  quite  accurately  estimate  the  length  and 
depth  alone.  Estimating  the  two  parameters  simultaneously  was  a  little  more  challenging;  however, 
we  still  obtained  reasonable  results.  In  estimating  length  alone,  using  either  500Hz  with  the  phase 
shifted  or  1kHz  with  raw  phase  produced  accurate  results  when  comparing  data  across  all  the 
data  points  across  the  center  of  the  damage.  In  a  few  cases  using  250Hz  or  2kHz  also  produced 
fairly  accurate;  however,  results  using  these  frequencies  were  not  consistent.  In  determining  the 
depth  of  the  damage  alone,  frequencies  500Hz,  1kHz  and  2kHz  again  produced  accurate  results 
for  the  depths  estimated  with  less  than  10%  error.  Our  tests  suggest  that  250Hz  is  not  a  viable 
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frequency  for  interrogation  in  any  case  presented  in  this  paper.  When  estimating  length  and  depth 
simultaneously,  using  500Hz  produced  accurate  results  in  many  cases.  Overall,  the  methods  again 
proved  to  be  feasible  when  used  with  experimental  data. 

In  some  cases,  depending  on  the  length  or  depth  of  the  damage,  certain  frequencies  were  found 
to  be  more  accurate  than  others;  therefore,  the  most  accurate  results  may  be  found  by  using  the 
reduced  order  POD  methodology  together  with  another  NDE  technique.  For  example,  the  magneto¬ 
optic/eddy  current  imager  ([16,  45])  is  a  visually  based  technique  which  displays  crack  images  or 
images  of  the  magnetic  fields  surrounding  the  actual  crack.  An  estimate  of  the  crack  length  may 
be  obtained  from  these  images  giving  us  a  priori  knowledge  of  the  damage.  Based  upon  these 
estimates,  a  source  frequency  as  well  as  an  initial  guess  for  the  optimization  routine  can  be  chosen, 
providing  valuable  information  for  use  in  the  inverse  problem. 

In  conclusion,  given  the  data  available,  we  have  shown  we  can  successfully  use  our  proposed 
reduced  order  computational  methodology  to  determine  both  the  length  and  depth  of  a  damage 
within  a  sample  (both  separately  and  simultaneously).  Moreover,  since  we  only  used  a  few  basis 
elements  (due  to  the  small  sample  size),  the  results  were  obtained  quite  quickly,  giving  us  a  method 
which  is  both  fast  and  accurate  on  experimental  data.  Indeed  this  section  gives  concrete  results 
indicating  that  the  POD  methodology  in  the  context  of  appropriate  least  squares  techniques  is  a 
practical  approach  to  nondestructive  damage  detection. 

7.  Concluding  Remarks 

In  this  paper,  we  developed  a  model  for  a  specific  eddy  current  method  making  some  simplifying 
assumptions  reducing  the  three-dimensional  problem  to  a  two-dimensional  problem.  We  utilized  a 
mathematical  tool  (phasors)  where  complex  valued  fields  were  employed  allowing  us  to  suppress  the 
time-dependence.  Furthermore,  for  computational  purposes,  we  included  two  additional  quantities 
in  the  Maxwell  formulations,  a  magnetic  vector  potential  and  a  scalar  electric  potential,  deriving  the 
boundary  value  problem  for  the  magnetic  vector  potential  with  some  additional  constraints  on  the 
electric  scalar  potential.  Given  the  magnetic  vector  potential,  we  could  easily  derive  the  magnetic 
flux  density  necessary  for  the  parameter  estimation  problem. 

We  then  presented  some  theoretical  results  which  established  the  existence  and  uniqueness  of 
solutions  as  well  as  continuous  dependence  of  the  solution  on  the  parameters  which  represent  the 
damage.  We  further  discussed  theoretical  issues  concerning  the  least  squares  parameter  estimation 
problem  used  to  identify  the  geometry  of  the  damage. 

Since  the  parameter  estimation  problem  involves  solving  the  forward  problem  numerous  times, 
we  needed  extremely  fast  and  accurate  solution  methods.  Therefore,  in  Section  4,  we  discussed  the 
reduced  order  POD  method  which  allows  one  to  create  a  set  of  basis  elements  spanning  a  data  set 
consisting  of  either  numerical  simulations  or  experimental  data.  The  POD  method  is  unique  in  that 
the  majority  of  information  is  captured  in  just  a  few  basis  elements,  allowing  us  to  use  a  smaller 
number  of  basis  elements  for  each  forward  solution.  This  results  in  a  substantial  decrease  in  total 
computational  time.  We  also  discussed  two  different  approaches  in  forming  the  POD  approximation, 
a  POD/Galerkin  technique  and  a  POD/Interpolation  technique  and  concluded  that  in  this  problem, 
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the  POD/Interpolation  method  gave  a  much  better  approximation  than  the  POD/Galerkin  method. 
Furthermore,  using  only  three  basis  elements  with  the  POD/Galerkin  method  we  still  obtained  an 
approximation  with  less  than  a  7%  relative  error  when  compared  to  a  finite  element  solution  using 
more  than  7000  finite  elements.  Using  the  POD/Interpolation  method,  we  achieved  less  than  a  1% 
relative  error  using  four  or  more  POD  basis  elements. 

In  Section  5,  we  presented  parameter  estimation  results  when  estimating  one  or  two  parameter 
values  using  simulated  data.  In  both  cases,  we  were  able  to  achieve  extremely  accurate  results  even 
with  10%  relative  noise  added.  In  addition,  on  average  we  obtained  a  total  reduction  in  time  of  a 
factor  of  approximately  103. 

Finally,  we  offered  results  of  the  parameter  estimation  problem  when  using  experimental  data 
obtained  from  a  giant  magnetoresistive  (GMR)  sensor.  The  experimental  results  were  based  on 
successfully  using  actual  experimental  data  to  form  the  POD  basis  elements  (instead  of  numerical 
simulations)  and  thus  illustrated  the  effectiveness  of  this  method  on  a  wide  range  of  applications. 
In  other  words,  whenever  it  is  difficult  to  model  the  physical  process  but  “good”  data  is  available, 
the  POD  method  may  be  a  viable  option.  Taken  as  a  whole,  our  work  here  indicates  that  using 
a  POD  computational  method  in  NDE  research  can  be  an  attractive  alternative  to  the  standard 
finite  element  methods,  offering  the  potential  for  substantial  savings  in  total  computational  time. 
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